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We investigate the ground state of the irrationally frustrated Josephson junction array with con- 
trolling anisotropy parameter A that is the ratio of the longitudinal Josephson coupling to the 
transverse one. We find that the ground state has one dimensional periodicity whose reciprocal 
lattice vector depends on A and is incommensurate with the substrate lattice. Approaching the 
isotropic point, A=l the so called hull function of the ground state exhibits analyticity breaking 
similar to the Aubry transition in the Frenkel-Kontorova model. We find a scaling law for the 
harmonic spectrum of the hull functions, which suggests the existence of a characteristic length 
scale diverging at the isotropic point. This critical behavior is directly connected to the jamming 
transition previously observed in the current-voltage characteristics by a numerical simulation. On 
top of the ground state there is a gapless, continuous band of metastable states, which exhibit the 
same critical behavior as the ground state. 



I. INTRODUCTION 

Frustration is regarded as a key concept in the under- 
standing of a number of complex cooperative phenomena 
in condensed matter systems including the glass transi- 
tion [l|, 0]. In general frustration leads to flattening of 
the energy landscape and suppresses the onset of conven- 
tional long-range orders, thereby opening possibilities of 
exotic phases. In order to understand the roles of frustra- 
tion, it is desirable to vary the strength of frustration in 
a systematic way. The frustrated Josephson junction ar- 
ray (J JA) under an external magnetic field [1] provides an 
ideally simple setting for this purpose since the strength 
of the frustration in the JJA can be tuned at will as dis- 
cussed below. The origin of frustration of this system is 
the impossibility to minimize the Josephson junction cou- 
pling energy of longitudinal and transverse bonds at the 
same time when a magnetic field is applied. This mech- 
anism seems to have a certain generahty. The frustrated 
JJA is closely related to the notion of frustrated crys- 
tals in curved space [l|, ■ Such frustrated crystals may 
be realized, to some extent, by actually bending crystals 
[1, H] . This amounts to injecting a given number density 
of dislocations into the crystals. Here the key parameter 
of the frustration is the number density of the externally 
induced dislocations. In the frustrated JJA, the disloca- 
tions correspond to the vortices induced by an external 
magnetic field. Quite remarkably, one can easily control 
the number density of vortices per plaquette / by just 
changing the strength of the applied external magnetic 
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field. 

By choosing an irrational / instead of a rational / Q 
we can maximize the strength of frustration in the sense 
that we can suppress the formation of a dislocation (vor- 
tex) lattice commensurate with respect to the underlying 
JJA. This system is called an irrationally frustrated JJA 
(IFJJA) and has attracted researchers for a long time. It 
has been observed numerically that the system remains 
in a vortex liquid state down to low temperatures 
and exhibits glassy signatures [l3, [lli • 

Recently it has been found that a relevant and very 
interesting perturbation on the IFJJA is the anisotropy 
of the Josephson coupling [l^ - [l^ . which breaks the bal- 
ance of longitudinal and transverse bonds and thus weak- 
ens the frustration. Anisotropic JJAs can be fabricated 
in laboratories by the lithography technique [l^. The 
IFJJA on the square lattice with different strengths of 
the Josephson couplings and Jy in the x and y di- 
rections, manifests itself as unusual vortex matter that 
slides freely in the direction of stronger coupling even 
at zero temperature, similarly to incommensurate charge 
density waves, but jammed into the direction of weaker 
coupling [T^]. It was argued that the mechanism of 
the sliding-jamming transition is similar to that of the 
Frenkel-Kontorova (FK) model [T6j . 

In the present paper, we study how the ground states 
(GSs) of the IFJJA change with the anisotropy A = 
Jy/Jx- To this end we develop an efficient numerical 
method valid even up to the isotropic point A = 1, which 
substantially extends the previous analysis based on a 
1 / A expansion [l3| , whose validity is limited to the strong 
anisotropy limit A 3> 1. The GSs are characterized by an 
incommensurate wave vector q(A) and its higher harmon- 
ics that continuously vary with A. With strong enough 
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anisotropy A ^ 1, the vortices are aligned as stripes par- 
allel to the weaker coupling axis [13|- By tuning the 
anisotropy smaller, the stripes become more tilted with 
respect to the weaker coupling (see Fig. [B We demon- 
strate that the so called hull function [1^, which re- 
flects the hidden periodic pattern that is incommensu- 
rate with the substrate lattice, becomes nonanalytic at 
the isotropic point A = 1. Physically this means that 
the sliding becomes prohibited in both the x and y direc- 
tions at the isotropic point, leading to a jamming tran- 
sition, which is analogous to the so-called Aubry tran- 
sition, which is well known in the FK model [1^1 ■ This 
analyticity-breaking transition is thus a static manifesta- 
tion of the observed jamming transition. 

To study the IFJJA with anisotropic Josephson cou- 
plings, we consider a classical model described by the 
following Hamiltonian with the Josephson coupling be- 
tween nearest neighbors along the x and y axes: 

i7(A,{</.J) = -^[cos0f -f Acos0f] (1) 

i 

with the gauge-invariant phase differences across the 
junctions, 

0, = (0^, 0f ) ^ - 0. - A^,e,+y -e.,- ) (2) 

Here Oi (i = 1, 2, . . . , N) is the phase of superconducting 
order parameter on the «-th island located at the lattice 
point Yi ~ {xi, yi) of a square lattice of size N — L x L 
with < a;^ < L, < j/j < L, Ti+i = {xi + l,yi) and 
Ti+y = (xj, + The parameter A denotes the strength 
of the coupling anisotropy. Due to the symmetry for 
permutation of axes x and y, it is sufficient to investigate 
only A > 1. The presence of the external magnetic field 
is described by the vector potential = {A^,A^). It is 
related to the filling factor /, which is the average number 
of flux quanta per square plaquettc, by taking its lattice 
rotation rot A, = 27r/ with rotX, = + Xf_^^ - Xf_^- - 

Specifically we employ a commonly used irrational 
number (3 - v^)/2 « 0.382 for / [I^. We believe, how- 
ever, that the properties we discuss below do not de- 
pend on the specific choice of the irrational number. To 
impose a periodic boundary condition, / must be ap- 
proximated by a rational number. It is known that a 
good approximation is obtained by using the Fibonacci 
series {F„ : K+2 = F„+i + F„} as / = F„_2/K. 
We use two series Fn = • • • , 8, 13, 21, 34, • • • and F,' = 
••• , 18, 29, 47, 76, • • • corresponding to different initial 
conditions (Fq = 1,Fi = 2) and {F^ = l,Fi = 3). We 
use system sizes L = F„ to ensure truly irrational filling 
in the thermodynamic limit L — > cx). 

II. VARIATION ANALYSIS 

We now study the GS of the system by extending the 
previous analysis that was limited to the computation 




FIG. 1. Vortex configurations of the metastable states (can- 
didates for the GSs) with various FRLV q = (g^ , /) in the 
anisotropic IFJJA with A > 1. White squares indicate the 
plaquettes containing vortices. The panels from left to right 
correspond, respectively, to q^L =13,14,15,16, and 17 with 
L = 34. As the anisotropy A increases, the GS changes from 
the left to right. Note that there is no recursive unit except 
the L X 1 strip and no simple relation between the tilt of 
apparent vortex stripes and the FRLV. 

of a few terms in the 1/A expansion [l3| . First we per- 
formed numerical searches of the GSs at various strengths 
of the anisotropy A on small systems using a simulated 
annealing method (not shown here). We found that the 
vortex configurations in the GSs exhibit stripe patterns 
whose tilt angle with respect to the axis of the weaker 
coupling varies with A (see Fig. [T|). This comes from 
the fact that </)j has one- dimensional periodicity, which 
is generally incommensurate with the underlying JJA, 
and a fundamental reciprocal lattice vector (FRLV) q(A) 
varies with A. We utilize this fact in the following anal- 
ysis. In the large anisotropy limit it was found to be 
qGs(oo) = (1/2, /) (shown in the right most panel in 
Fig. [ID El. 

Let us now explain our strategy, which is slightly re- 
formulated but equivalent to the scheme presented in 
Ref. [H except for some partial use of numerical pro- 
cedures. First, recalling Eq. we find that the phase 
differences must satisfy the following equation at each 
plaquette: 

rot^, = -rotAi = -27r/. (3) 

To meet this condition, it is convenient to decompose </> 
into a uniform-rotation part and a rotation-free part as 

"max-l 

ct>, = <P:+ £ c^„e^"-^'i--. (4) 

where 

rot0* = -2^/ rot [^^e^^^'^^''^'] = 0. (5) 

The upper bound of the harmonics rimax is given by the 
smallest n with which both nq^ and ng^ are integers. It 
almost always equals L in the present model and becomes 
infinite in the irrational limit. The first of Eqs. ([5]) can be 
solved by (f)* = {—2iTfyi + Cq, 0), where Cq is a constant 
introduced to satisfy the periodic boundary condition for 
9i. The second of Eqs. ([5]) can be solved by </3„ = ((/j^, (pj^) 
with ^l/^l = (1 - e2"''*9'')/(l - e^"'^*'?"), which can be 
used to eliminate the degrees of freedom {(/S^} in favor of 
{•^n} (o'^ ■^ice versa). 

We assume that the FRLV can be parametrized as 
Q = (5^1/) by a parameter g^', which is suggested by 
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the 1/A expansion analysis [l^l- This perturbative anal- 
ysis suggests that a unique stable solution, which satisfies 
the current conservation (force balance) condition, can be 
constructed explicitly for a given that is left as a con- 
trol parameter. The analytic computation is, however, 
quite cumbersome and difficult to continue to higher or- 
ders. To overcome the difficulty, we numerically search 
the configuration with minimal energy, 

£'m(A,q) = miniJ(A,q, {v„}), (6) 

for a given FRLV. We performed this analysis in the 
range f < < 1/2. 

For numerical energy minimization, we employ the 
steepest-descent method for the degrees of freedom {tpn}. 
Actually, we numerically integrate the equations of mo- 
tion; dipfjdt = --dH/d(p^ for 1 < n < nmax, with an 
initial condition 1^9^ = for all n. We determine that 
the system relaxes to the minimum-energy state when 
\dH/dt\ becomes less than 10~^. As an example we show 
in Fig. [T] the vortex configurations corresponding to the 
candidate FRLVs for L = 34. 

The final step is to determine qGs(A): q which gives 
the minimal energy for given anisotropy A as 

£;gs(A) -minS„(A,q) =£;™(A,qGs(A)). (7) 
q 

Although there is no rigorous proof that the states ob- 
tained by this method are the true ground states, we find 
that they coincide with the ones obtained by the simu- 
lated annealing methods. Wc thus believe that they are 
the true ground states. 
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III. RESULTS 

A. successive structure transition 

Figure [2ja) shows the minimal energy i?,„(A,q) as a 
function of q^ for A=1.0, 1.5, and 2.0. From this we 
determine qos and Eqs for a given anisotropy A. The 
energy behaves quadratically in the vicinity of each mini- 
mum as £',„ ( A, q) - Eos (A) = L^c( A) [q- qcs ( A)] ^ where 
c(A) is a constant of 0{L^). Thus there is a continu- 
ous band of mctastable states with slightly different FR- 
LVs around the GS. (We confirmed that the obtained 
metastable states are stable even if we lift the constraint; 
the Fourier components equal zero except for the har- 
monics nq.) 

Figure [^b) shows the A dependence of the x compo- 
nent of the qGs(A) = {QgS' f) with a ^ x,y. It can be 
seen that q^g monotonically increases with A and changes 
in a stepwise manner at several points, where level cross- 
ings occur between metastable states with neighboring 
FRLVs. As L becomes larger, the number of steps in- 
creases and qQg tends to be a continuous function of A. 
It can be seen that the horizontal stripe state q = (1/2, /) 
is no longer the GS for A < 1.8. The vortex stripes 



FIG. 2. (Color online) (a) Minimal energy vs for various 
values of A. Constant values are added in the cases of A = 1.0 
and 1.5 (by 0.155 and 0.050, respectively) for a convenient 
view, (b) Anisotropy dependence of the FRLV of the GS, 
qcs- 



in the GS become more tilted with respect to the weaker 
coupling axis when approaching the isotropic point A = 1 
(see Fig. [1]). 

The GS of the isotropic system (A = 1) is of particular 
interest. For small sizes L < 18, the GS is the previ- 
ously found staircase state [l7l - [l9| with ggg = q^g = f- 
However, this is not the case for L > 34. A new min- 
imum appears at q^ « 0.416 and there is a local min- 
imum at q^ w 0.388, both of which are larger than 
/ = 0.381966 • • • . We obtain g5g=14/34, 23/55, 37/89, 
60/144, and 97/233, where the numerators constitute a 
Fibonacci series (with Fq = 1 and F" — 4) and the 
denominators correspond to L. Therefore, q^g seems to 
equal 0.4164078 • • • [> (3- V5)/2] in the limit of L ^ 00. 
A notable feature is that the symmetry for the permuta- 
tion of X and y axes spontaneously breaks in contrast to 
the symmetric staircase state mentioned above. 
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FIG. 3. (Color online) Hull function of (a) 0^ and (b) c^" for 
six different values of A. Each graph has L(=144) points. The 
arrows indicate the direction of decreasing A. 



B. breaking of analyticity in hull functions 

Next we investigate more closely the properties of the 
GSs parametrized by the FRLVs q, which are generally 
incommensurate with respect to the underlying JJA. A 
useful measure of such an incommensurate object is the 
so called hull function [IlHIlOl 4>{z) = (<^^(z), 
with which the phase field is written as 

= 0(qTO. (8) 

In Fig. [3] the hull functions (p^{z) and (j)^{z) are plotted 
with respect to < z < 1. It is noteworthy that (j)^ is 
distributed, covering all phase from — tt to tt, while (j)^ 
is bounded around zero to have a gap around i/)^ — ±7r. 
When the hull function is smooth and gapless, there is a 
sliding soft mode [l^ ; the infinitesimal increment of 0^ 
is compensated for by the shift in the argument of the 
hull function q • r,;, which costs no energy barrier. Since 
this leads to dissipativc vortex flow driven by current 
induction, the GSs have superconductivity only along the 
stronger-coupling direction. 



I ■ I ■ I ' I ' I ' I 
0.0 0.2 0.4 0.6 0.8 1.0 

z 

FIG. 4. (Color online) Derivative of the hull function of 0". 
(a) Anisotropy A dependence for L = 144. (b) Size L depen- 
dence at the symmetric point A = 1. 



A remarkable feature is that the hull function in both 
the X and y directions becomes distorted by decreasing 
the anisotropy A. While the hull functions are smooth 
with larger anisotropies [isj . steplike structures appear 
when approaching the symmetric point A = 1. The 
hull function at A = 1 appears to be nonanalytic: It 
is discontinuous at several points. Figures Hl^a) and 
Sljb) show the discrete derivative of the hull functions 
d4>"{z)/dz = [0"(z + 2tt/L) - ^"(z)]/(27r/L). It can be 
seen that sharp peaks emerge when approaching A = 1. 
This is reminiscent of the critical behavior of the A ubry 
transition in the Frenkel-Kontorova model EB, l20l |. 
Figure Hfb) shows the size dependence of the derivative 
at A = 1. As the size L increases, not only do the peaks 
become sharper, but also the number of peaks increases 
reflecting higher harmonics. We speculate that the hull 
function becomes discontinuous everywhere in the ther- 
modynamic limit. 
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C. Scaling behavior and diverging characteristic 
length 

Let us examine further the singular behavior around 
the symmetric point A = 1. Figure [5l[a) shows the ampli- 
tude of the nth harmonic nq [note that q is set equal to 
qGs('^ = 1) independently of A]. At A = 1, |(^^| seems 
to decay as a power function of n with exponent close to 
— 1 although there are many dips and peaks. In contrast, 
higher harmonics exponentially decrease above a certain 
characteristic scale n*{\) for A > 1. This cutoff indicates 
the smoothing of the hull functions. We found a scaling 
behavior 



n*(A) 



with 7i*(A) cx |A - 1]-'' (9) 



not only for A > 1, but also for A < 1, as shown in 
Fig. [5{b). The exponent v is roughly estimated as ~ 
2.4. We confirmed that the same scaling works well for 
the FK model with v = 0.9874 [HI, [22[. This scaling 
means that a diverging number of higher harmonics come 
into play when approaching the symmetric point A = 1. 
We note that this observed singularity is present not only 
for qGs(l), but also for other q's around qGs(l): i-e., for 
the continuous band of metastable states. 

Initially, it would appear as if the characteristic length 
in Eq. ([9|) were l/n*|q| and it went to zero, which is con- 
trary to ordinary critical behavior. However, we specu- 
late, inversely, that there is a diverging length scale for 
the following reasons. On lattice systems, the compo- 
nent of the wave vector, nq, should be treated in the first 
Brillouin zone (-1/2, 1/2) as g™^ = nq - \nq + 1/2J, 
where [• • • J is the floor function. For irrational g, {g™^} 
behaves like a series of random numbers made by the lin- 
ear congruential generator in a long span of n. Thus the 
maximum of 1 /(?™^ for n < 7i* is roughly proportional to 
n* . If g^ and q^ are independent irrational numbers, the 
maximum of the vector length is proportional to -v/n*. 
In addition, the number of harmonics in the finite-size 
system is bounded by L when g is approximated by an 
irreducible fraction p/L (note that g™| = g™^)- This 
also indicates that n is a measure of the length. 

The breakdown of the analyticity of the hull function 
means that the sliding soft mode disappears; jamming 
occurs even in the direction where the hull function is 
gaplcss. Indeed, of the transport properties of the present 
system it has been found that the current-voltage charac- 
teristics exhibit a scaling feature around the symmetric 
point A = 1 [l3l. Thus Eq. ^ is regarded as a static 
signature of the jamming transition. 
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FIG. 5. (Color online) Spectrum of the harmonics of the hull 
functions (a) before and (b) after scaling. We set u = 2.4. In 
the scaling, we do not use the data for n > I//4, which show 
a strong finite-size effect. 



with the variation of the anisotropy A of the Joseph- 
son coupling between the horizontal vortex stripes in the 
strong anisotropy limit and nearly (but not exactly) di- 
agonal vortex stripes at the isotropic point A = 1. 

It is interesting to note again that the present system 
admits a continuous band of metastable states around 
the GS, which is presumably the source of the glassiness 
in the present system. The analyticity breaking of the 
hull function approaching the isotropic point A — > 1 im- 
plies jamming of the sliding vortices. In addition, the 
scaling behavior of the harmonic spectrum suggests the 
existence of diverging static length scale. It would be 
very interesting to study the properties of the present 
system at finite temperatures from static and dynamic 
view-points. 



IV. CONCLUSION 

To summarize, we studied GSs and low lying 
metastable states of the IFJJA with anisotropic Joseph- 
son couplings. We found the GS changes continuously 
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